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Summary 


A method  is  presented  for  the  estimation  of  the  param- 
eters in  the  vector  autoregressive  moving  average  time  series 
model.  The  estimation  procedure  is  derived  from  the  maximum 
likelihood  approach  and  is  based  on  Newton-Raphson  techniques 
applied  to  the  likelihood  equations.  The  resulting  two-step 
Newton-Raphson  procedure  is  computationally  simple,  involving 
only  generalized  least  squares  estimation  in  the  second  step. 
This  Newton-Raphson  estimator  is  shown  to  be  asymptotically 
efficient  and  to  possess  a limiting  multivariate  normal 
distribution. 


Key  Words:  Vector  autoregressive  moving  average  process,  Newton- 

Raphson  method,  maximum  likelihood  estimation. 
Generalized  least  squares. 
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1.  INTRODUCTION 


In  the  past  several  years,  following  the  work  of  Box 
and  Jenkins  (1970)  and  Hannan  (1969b, 1970) , there  has  been  a great 
deal  of  interest  in  the  use  and  estimation  of  autoregressive  mov- 
ing average  models  (ARMA  models)  for  a scalar  time  series  of  obser- 
vations. The  vector  ARMA  model  has  received  much  less  attention. 
Hannan  (1969b,  1970)  has  considered  the  estimation  of  the  pure 
moving  average  model  in  the  vector  case  using  spectral  methods. 
Nicholls  (1976a)  has  extended  this  work  to  the  case  of  estimation 
of  the  mixed  vector  autoregressive  moving  average  model  which  con- 
tains exogenous  variables,  the  so  called  VARMAX  model.  An  iter- 
ative nonlinear  least  squares  (maximum  likelihood)  estimation  pro- 
cedure for  the  stationary  vector  ARMA  model  has  been  presented  by 
Tunnicliffe  Wilson  (1973).  Others  who  have  considered  the  esti- 
mation of  vector  ARMA  time  series  models  include  Akaike  (1973), 
Kashyap  (1970),  and  Whittle  (1963).  However,  to  a large  extent, 
the  problem  of  estimation  of  mixed  vector  models  has  been  con- 
sidered too  computationally  complicated  to  be  practically  useful. 
Also,  the  estimation  theory  is  much  more  difficult  than  in  the 
scalar  case.  (See  the  recent  paper  by  Dunsmuir  and  Hannan  (1976) 
for  a treatment  of  the  estimation  theory  of  more  general  vector 
linear  processes  which  include  the  ARMA  models). 

In  this  paper  a method  of  estimation  for  vector  ARMA 
models  will  be  presented  which  is  asymptotically  efficient  yet 
computationally  simple.  The  method  to  be  proposed  uses  the  maxi- 
mum likelihood  approach 


and  is  based  on  Newton-Raphson  techniques  applied  to  the  likelihood 
equations  (in  the  time  domain).  The  resulting  "Newton-Raphson" 
estimator  is  shown  to  be  asymptotically  equivalent  to  the  maximum 
likelihood  estimator  and  to  possess  a limiting  multivariate  normal 
di stribution. 

Thus  we  consider  the  estimation  of  the  parameters  in  the  model 

Yt  " ^ AiYt-i  = 6t  + iJ1  Gi€t-i’ 

where  Yt  = (ytl,. • • >Vtm) ' and  ^ = ( 6 tl , • • • , e tm) ' are 
of  m components,  and  A^,...,Ap,  G^,...,G  are  (m * m) 
the  unknown  parameters.  With  the  use  of  the  lag  operator 
that  = Yfc  this  vector  model  may  be  written  as 

A'(^)Yt  = G'K)€t, 

where  A(/)  = I - A±£  - ...  - A pi£P 

and  GU)  = I + G±id+ . . . + Gq^q  . 

The  are  assumed  to  be  independent  and  identically  distributed 

with  mean  vector  0,  covariance  matrix  I,  and  finite  fourth 
order  moments.  We  also  assume  that  all  roots  of 

det  A(z)  =0  and  det  G(z)  =0 

are  greater  than  1 in  absolute  value  and  there  are  no  common  roots 
to  the  two  equations.  This  assumption  implies  that  the  matrix  series 

A(z)  x = Z J.z1  and  G(z)"x  = Z D,  z1 
i=0  1 i=0  1 

both  converge  for  |z|  < 1 + a,  A > 0 , and  that  Yt  can  be 


(1.1) 

vectors 
matrices  of 
such 

(1.2) 


2 


expressed  as  the  stationary  solution  to  the  difference  equation 
(1.1)  as 


Yt  = a7  urV  U)  ^ 


oo  q 
2 Z 
u=0  v=0 


W^t-u-v* 


The  problem  of  identification  for  models  of  the  form  (1.1)  has 
been  discussed  in  Hannan  (1969a,  1975),  and  we  assume  that  the 
identifiability  conditions  given  there  are  satisfied.  In  particular, 
we  assume  that  -the  rank  of  (A^  I G^)  is  equal  to  m. 


2.  THE  METHOD  OF  ESTIMATION 

For  the  estimation  of  the  parameters  in  the  model  (1.2), 
we  assume  that  we  have  the  vectors  of  observations  , t=l,...,T. 
We  define  the  (T-p)  xm  matrices 


, ^iY- 

( w) 

, 6- 

! €p+i  ' 

K+i-\ 
| . 

/ 

1 

Is*  / 

\€/T-i  / 

Then  we  can  write  the  entire  model  in  matrix  form  as 

P i 9 i 

Y - 2 tC~YA.  = <E  + Z ^60..  (2.1) 

i=l  1 i=l  1 

Alternatively,  we  may  express  the  model  more  conveniently  in  vector 
form  using  the  theory  of  Knonecker  products  of  matrices.  For  any 
(mxn)  matrix  C,  let  vec(C)  denote  the  (mn  xl)  column  vector 


- 3 - 


formed  from  C by  placing  the  columns  of  C below  each  other 
starting  from  left  to  right.  Then  a result  which  we  shall  use 
to  write  the  model  in  vector  form  is  that 

vec(ABC)  = (C'  ® A)  vec(B), 

where  A,  B,  and  C are  any  matrices  such  that  the  matrix  product 
ABC  is  defined.  (See  Neudecker  (1969)  for  a proof  of  this  and 
other  results  for  Kronecker  products.)  Defining  the  vectors 
y = vec  ( Y) , e=vec(€),  ^y  = vec  (a^Y) , s^e  = vec  (^i€) , ai=vec(Ai), 
i=l,  •••,?,  and  yi=vec(Gi),  i»l,***,q,  we  can  write  the  entire 
model  in  vector  form  as 


Pi  Q i , 

y - Z (I_®S£1Y)  a,  = e + I (I  ®^€)y±  • (2.2) 

i=l  ra  1 i=l  m 1 

To  motivate  the  estimation  procedure  we  suppose  the 
are  normally  distributed  and  use  the  maximum  likelihood  approach.  To 
simplify  the  form  of  the  likelihood  function  certain  assumptions  will 
be  made  concerning  the  initial  observations  and  disturbances.  First, 
we  consider  the  initial  observations  , . . . , Yp  as  fixed,  and  estimate 
from  the  likelihood  function  conditional  on  these  values.  Second,  we 
assume  that  the  initial  disturbances  c1_  are  equal  to 

their  unconditional  expectations,  which  are  0.  Then  introducing  the 
(T-p)  x (T-p)  lag  matrix  L which  has  l's  on  the  diagonal  directly 
below  the  main  diagonal  and  0's  elsewhere,  we  define  the  (m(T-p)) 
x (m(T-p))  matrix 


^ = Pm  1 


ra(T-p)  + £ ■ 
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-1  T-P'1,  i 

with  det  G = l.  (Note  that  G = Z (D'®L  )).  Thus  we  can 

i=0  1 

express  the  entire  (modified)  model  in  vector  notation  as 


y - Z (I  ®*iY)a1  = 
i=l  m 1 


q i 

e + Z (I  ® Lfc) 


i=l 


m 


Yi 


q i 

e + Z ( G ' ® L ) e 
i=l  1 


Ge 


(20) 


Based  on  the  assumption  of  normality  of  the  q t,  the 
(modified)  likelihood  function  of  the  observations  Yp+i » • • • • 
conditional  on  Y^,...,Y  , is 


F = 


(27r)*mN|z(£N 


exp 


/ 

G'-1 


• (Z-1®IN)G-1{y-Ji(Im®^iY)a1j]  , 

where  N = T-p.  Thus  we  obtain  the  partial  derivatives 

- (Im®^Y)'0'-1(Z-1®IN)G-1[y-ili(Im®/Y)al]  , 


( J — 1 ) • • • jP)  ) 


SlogF 
3Y 


Si—  = (lm®LJ€)/G/'1(Z‘1®IN)G-1  |y-Z  (lir)®sdiY)a1}  , 
J i~l 


( «J  — 1 > • • • t q ) > 
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where  € is  expressible  in  terms  of  the  observable  quantities 
Y^  through  equation  (2. 3)  as 

e = vec(€)  = G"1(y-^(Im®£iY)a1^  . 

Defining  the  vector 


9 — ( Ct±  y . 


= vec(A1| . . . |Ap|G1| . . .Gq)' 
and  the  matrix 


W = ((Im®=£Y),...,(Im®£PY),(Im®L€),...,(Ira®Lq6)), 
we  can  express  these  derivatives  collectively  as 

= W/G/-1(2-1®I  jG^fy-l  (I  ® £iY)ai  !•  . (2.4) 

90  i=l  m 

The  maximum  likelihood  equations  which  result  when  (2.4)  is  set 
equal  to  0 are  nonlinear  in  the  parameters  9 (unless  q = 0, 
in  which  case  the  solution  is  just  the  least  squares  estimate 
§=(w'w)  W'y  which  is  nearly  identical  to  the  Yule-Walker 
estimate  obtained  in  (a)  below).  Thus  for  q>0  these  equations 
can  only  be  solved  by  numerical  procedures  such  as  the  Newton-Raphson 
method.  The  Newton-Raphson  equations  for  an  approximate  maximum 
likelihood  estimator  9 are 


92logF 


9099' 


(8-en)  = 


99 


(2.5) 
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where  0Q  is  an  initial  estimate  of  0.  It  can  be  shown  that 
on  neglecting  terms  which,  when  divided  by  T,  converge  to  0 
in  probability  as  T -*■<»,  we  will  have  the  approximation 

- = W/G,~1(Z~1»Iw)G'1W  . (2.6) 

d9de' 

To  obtain  the  Newton-Raphson  estimator  of  0,  we  assume 
that  we  have  initial  estimates,  0Q  = . . . , y^) ' and  % , 

which  are  consistent  estimates  of  9 and  Z to  the  order  T-^  in 
probability.  These  estimates  may  be  obtained  as  follows: 

(a)  Compute  the  sample  autocovariances 

T-s 

C( s)  = — Z (Yt  - Y)(Yt+s  - Y)'  = C'(-s),  (s=0,l,...,p+q), 

T 

1 T 

where  Y = — Z Y. . We  then  estimate  the  parameters  a.  = vec(A, ) 

T t=l  L 1 

consistently  by  solving  the  vector  Yule-Walker  type  equations 

P 

C(-s)  - Z £(c(i-s)  = 0,  ( s=q+l, . . . ,q+p) . 

i=l  1 

(b)  Having  obtained  the  estimates  ~A. , we  form 

P P 

C ( s)  = Z Z A'  C(s+J-X)A.  = C ( -s) , (s=0,l,...,q), 

j=0  x=o  J 

where  AQ  = -Im,  and 

f„(X)  = — S C (s)e"isx,  -7 r<  X < 7 r. 

2tt  s=-q 

These  are  consistent  estimates  of  the  autocovariances  and  spectral 
density  matrix,  respectively,  of  the  moving  average  process 
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V.  = € . + Z G(  €+  . . If  f ( X ) is  an  Hermitian  nonnegative 
definite  matrix,  i.e.,  fv(\)  > 0 for  -tt<X<tt,  it  may  be 

factorized  in  the  form 


?v(\)=—  G' (eU)Z  S (e~U), 

2t r 

to  obtain  consistent  estimates  = vec(cL),  i=l,...,q,  and  IT. 
Robinson  (1967,  pp.  190-200)  describes  a computer  program  for 
such  a factorization.  Consistent  estimates  of  the  G^  may  also 
be  obtained  by  a method  which  does  not  require  the  factorizing  of 
the  spectral  density  matrix  ?v(X)  (see  Hannan  (1975),  page  977, 
equation  (5))* 

Now  we  let 


Then  substituting  (2.4)  and  (2.6)  into  equation  (2.5),  we  obtain 
the  following  Newton-Raphson  equations  for  §, 

W/G/'1(Z“1®  In)G_1W(0  - 0Q) 

= w'cr-^r1®  in)g_1  fy  ■ l2i(im®^lY)aij  • 

P , i _ 

Since  fl  0O  + y - Z ( Iffl  ® Z1Y)aj_ 

- y + > 
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the  Newton-Raphson  equations  can  be  written  in  the  form 


ft' S'  -1(2'“13In)G~1W9 


The  Newton-Raphson  solution  9 to  equation  (2.7)  can  be  interpreted 
as  the  generalized  least  squares  solution  to  the  identity 

y + 2 (lm®Li€)y 
i=l  m 1 

p , 1 q 

= 2 (l®^lY)a.  + 2 (I  ®L1€)y1  + Ge 

i=l  m 1 i=l  ra  1 


+ 


(2.3) 


where  the  last  term  on  the  right  hand  side  of  the  equation  is  to  be 
neglected  and  the  error  term  Ge  is  treated  as  having  covariance 
matrix  G(2  3 I ^ . 

The  asymptotic  properties  of  the  estimator  9 will  be  given 
in  the  next  section.  Before  proceeding  to  these,  we  conclude  this 
section  by  briefly  discussing  the  computations  that  are  needeu  to 
complete  the  estimation  procedure  and  obtain  § as  the  solution  to 
equation  (2.7).  Once  the  initial  estimate  9^  has  been  obtained, 
the  vector 

- if  P i 

e = vec(?)  = (T  |y  - 2 (Im®^  Yjc^ 

P , i __ 

may  be  computed  recursively  from  G e = y - 2 (I  3^  Y)a,  as 

i=l  m 1 

”.2/,AiYt-i  “ 2Gi  ^t-i’  ( t=p+l,  . . . , T) , 

1=JL  1=1 
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where  € „ =...=€  =0.  Similar  calculations  may  be  used 

p+l-q  p 

to  obtain  the  columns  of  the  matrix  of  "independent"  variables 

W = G~^W  recursively  from  GW  = W, 
and  the  vector  of  the  "dependent"  variable 

? = rlfJ+  j1(l-®Ll?)7i|  recursively  from 

~ - <1  i- 

G y = y + Z (I_®  LX€)y.  . 

i=l  m 1 

Then  the  estimator  9 is  simply  the  generalized  least  squares 
solution  to  the  regression  of  y on  W,  i.e., 

^ « { W'  (2"1  ® I^w}"1  W' (Z_1®IN)y  . 

5.  THE  ASYMPTOTIC  DISTRIBUTION  OF  THE  ESTIMATOR 

To  describe  the  asymptotic  distribution  of  the  estimator  9 
we  introduce  the  matrices  Y,$,  and  Cl,  whose  (g,h)th  submatrices 
will  be  defined  below.  In  expressing  the  following  integrals,  for 
notational  convenience  we  will  omit  the  argument  variable  e1^  from 
the  matrix  polynomials  A(e^)  and  G(e^)  defined  in  Section  1, 
so  that  these  will  appear  below  simply  as  A and  G.  Thus  we  define 
the  submatrices 

= J ei^g"h)XG"1Z"1G*_1 ® ( A*-1G*ZGA”1) 7 dX, 

27 T -7T 

(g»h=l, • • • ,p). 


• * 


m 


10 


Vh  * I E {(lm®Ls€)'G'-1(Z-1®IN)0-1(lm®Lh€)| 


1 

f 

ei(g-h)XG-l2-l 

2t r 

-7T 

lim 

±e 

|r(lm®^Y)'G/ 

T-*a> 

T 

1 

r 

ei(g-h)XG-lz-! 

2tt 

•t r 

( S— 1 1 • ■ • >P5  h— 1 * • • • > q)  • 


Then  we  have  the  following 


THEOREM.  For  the  model  (l.l)  under  the  assumptions  stated  in  Section  1, 
let  9 denote  the  estimator  of  0 = vec(A^ | ••• | A | G^ !••• IG^) ' as 
obtained  from  equation  (2.7)*  Then  the  distribution  of  -\/T(9-9) 
converges  to  a multivariate  normal  distribution  as  T+oo,  with  mean 
vector  0 and  covariance  matrix  equal  to 


Y a 

O'  i 


(5-1) 


Before  considering  the  proof  of  the  theorem,  we  make  the 
following  comments: 

(1)  The  matrix  V defined  in  (5-1)  is  equal  to 


11 


lim 

T->oo 


lira 

T-*« 


^(z-1®  In)G_1w|  , 


_ 1 E(  AogF)  _ 

t aeae' 


the  "asymptotic  information"  matrix  of  0.  Thus  the  asymptotic 
distribution  of  § is  identical  to  that  of  the  maximum  likelihood 
estimator  based  on  the  assumption  of  normality  of  the  disturbances 
so  that  9 is  asymptotically  efficient  in  that  sense.  Also,  9 
converges  to  9 in  probability  as  T-^oo. 

(2)  If  we  let  denote  the  (i,J)th  element  of  Z-1,  then  we 


have 


/ 


dlogF  ) 

17?  - 


| oii  - | e' (Eii ® IN)e  for  i=j, 

N °ij  “ \ e'  {(EiJ  + EJi)  ® e for 


where  e = G~  | y - Z ( Im  ® £iY)ai|  and  E^.  denotes  the  (mxm) 

matrix  with  1 in  the  (i,j)th  position  and  zeroes  elsewhere.  We 
have  used  the  fact  that  for  a nonsingular  matrix  Z, 

for  i=j 

for  i^j  . (see  Golderberger  (1964), 

Chapter  1) 

Thus  an  asymptotically  efficient  estimator  of  Z is  given  by 

2 = 1 rt  , 

N 

where  e = vec(£)  denotes  the  vector  e above  with  estimates  9 


doiJ 


'ii 


2o 


ij 
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used  in  place  of  9.  This  estimator  1 is  asymptotically  uncorrelated 
with  Similarly,  the  covariance  matrix  of  § may  be  estimated  by 

■[&'&'  "1(i“1  ® iN)&_1w}  ~1  . 

(3)  The  estimation  procedure  for  5 can  easily  be  modified  to 
include  the  case  where  elements  of  the  autoregressive  moving  average 
coefficient  matrices  A^,...,A  , G„ , . . . , are  specified  to  be  zero. 

We  simply  modify  the  definitions  of  the  matrix  G and  the  vector  e 
according  to  these  restrictions  (although  this  is  not  necessary  for 
the  results  of  the  theorem  to  hold),  and  in  the  matrix  W we  delete 
the  columns  of  (lm®£iY),  1=1,..., p,  and  (lm®L'!'?),  i=l,...,q, 
corresponding  to  the  zero  elements  of  = vec(Ai)  and  yi  = vec(G^), 
respectively.  The  resulting  estimator  obtained  from  equation  (2.7) 
will  again  be  asymptotically  efficient,  consistent,  and  will  possess 
an  asymptotic  normal  distribution  with  covariance  matrix  similar  to 
V (but  with  the  appropriate  rows  and  columns  of  V deleted). 

PROOF  OF  THEOREM:  The  proof  of  the  theorem  is  much  the  same  as  that 

of  a similar  theorem  for  the  scalar  case  given  in  Reinsel  (1976a), 
and  so  we  shall  omit  most  of  the  details.  We  ignore  the  effect  of 
the  modification  of  the  initial  disturbances  €t,  and  hence  the  use 
of  the  matrix  G in  equation  (2.8),  since  the  modification  has  a 
negligible  effect  as  T+oo  and  does  not  affect  the  asymptotic 
. properties  of  the  estimator.  Then  using  the  (modified)  identity  (2.8) 
and  equation  (2.7) t we  have 


I V?  (3-0)  = W'G,_1(2:_1®  IN)2T1Wj  • -^-  W'3"1^*1®  IN)e 

. 

+ |i  ’X  (5.2) 

• J ^ ff/S/"1(ST1»  IN)^_1  (lra®  (L^-L^jJ'VT^-Yi). 

Now  it  can  be  shown  that 

plim  - W'g'  “1(2_1  ® Im)G-1W  = lim  m E (w'  G'  '1  (Z'1  ® IM) G_1W>  = V. 
T-»oo  T N T->oo  1 1.  N ‘ 

Also,  as  in  Reinsel(1976a) , it  can  be  shown  that  each  of  the  summands 
in  the  second  term  on  the  right  side  of  (5*2)  has  a probability  limit 
equal  to  0 by  the  consistency  of  the  initial  estimates  (and  the 
boundedness  in  probability  of  the  vectors  VT  (y^-Y^))*  And 
finally,  again  by  the  consistency  of  the  initial  estimates,  from  (5-2) 
it  is  not  hard  to  see  that  the  asymptotic  distribution  of  -/T  (9-9) 
will  be  the  same  as  that  of  V-1  * — w' G' ~1(Z~1 ® IN)e  . 

v; 

Thus  it  follows  that  the  results  of  the  theorem  will  be  established 

once  it  is  shown  that  the  vector  W'  G'  (Z  ®IM)e  has  a limiting 

-y/T  N 

normal  distribution  as  T -*•  « with  mean  vector  0 and  covariance 
matrix  equal  to  V. 

To  establish  the  asymptotic  normality  of  this  vector,  we 
consider  in  detail  the  asymptotic  behavior  of  a subvector  of 


i.  / / — i.  — i 

W' G'  (2  For  example,  we  have 

VT 


1 T-p-1  /■  * 11/1 

= 2 iD.2‘i8(m)  f 

VT  j=0  (•  J J 


= I vec  /(iA^Y) ' € 2*1D'  ] 

VF  j=0  <■  J J 


1 T-p-1  T 

: — 2 2 ve^Y*, 

V?  j=0  t=p+j+l  Z 1 3 t J 


1 T t-p-1 

— 2 2 vec(Yt  , Vt2  Dl) 

VT  t=p+l  j=o  t-i-j  t J 


1 T t-p-1  oo  q 

: 2 2 2 2 vec( J G €.  . , V2  D') 

VT  t=p+l  j=0  u=0  v=0  u V t-i-J~u-v  t j/ 


1 T fmin( t-p-1, n)  n q 
: — — 2 < 2 2 2 vec( 

VT  t-p+1 1 j=0  u=0  v=0 


(5-3) 


ZTn  + ^Tn* 


for  n < T-p-1, 


where 


T 

Z 


W._, 


JTn  “vt"  t=;+i  tn 


(3.^) 


w 


tn 


min(t-p-l,n)  n q 1 

2 ulo  v=0  VeC(JuG;et-l-j-u-vetZ"  Dj)-  (3-5) 


J=0 


and 


( t— p+1 


. T t-p-1  oo  q 1 , 

*m  *4-  , s E , * z v«(Wt-i -j-u-ve;2'  V 

Vt  t=p+n+2  J=n+1  u=0  v=0  0 J 


1 T min(  t-p-1, n)  oo  q 

— Z Z Z Z 

VT  t=p+l  j=0  u=n+l  v=0 


(3.6) 


Now  since  the  elements  of  the  matrices  Ju  and  Du  converge  to  0 
exponentially  as  u + oo,  it  is  not  hard  to  show  that  the  remainder 
term  converges  to  0 in  probability  as  n-*oo,  uniformly  in  T 

(e.g.,  it  can  be  shown  that  E(Rix^R<rn)  0 as  n->oo,  uniformly  in  T). 
For  fixed  n,  the  vectors  Wtn  have  0 means,  are  uncorrelated,  and 
the  covariance  matrix  of  Wt  is  equal  to 

E(«tnMto>  ‘ vtn 

min( t-p-1, n)  n q 

= Z 

J, k=0  u,m=0  v, s=0 


S I SB  (vecUXSt-i-J-u-v^"1^) 


vec 
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min( t-p-l,n) 
2 

j,k=0 

n 

2 

u,m=0 

q 

2 

v,  s=0 

E JuGv>€t-t-J-u-veU-k-m- 

min( t-p-l,n) 
2 

j,k=0 

n 

2 

u,  m=0 

o 

II 

w 

> 

E { (DjZ'X  * JX)*(V'S  ® Jm°s>  'j 

min( t-p-l,n) 
2 

j,k=0 

n 

2 

u,  m=0 

q 

2 

v,  s=0 

E{(Dj2'l6tei2'X®JuG^°sJm)} 

min( t-p-l,n) 
2 

j,k=0 

n 

2 

u,  m=0 

<5 

to  M.O 
II 

o 

('A'W.V'  (3-7> 

s 


where  the  summation  extends  over  only  those  indices  for  which 
j+u+v  = k+m+s.  Thus  the  covariance  matrix  of  ZTn  is 


E(zTnzii>  ’ - 2 

in  m T t=p+l 


Vtn'  and 


limi  E(zTnzTn) 
1 00 


= lim  - 2 


min(t-p-l,n)  n q 1 

Z £ £ (DZ-X®^Vm) 

vm r\  r o — r\  w 


T-*»  T t=p+l  j,  k=0  u,m=0  v,  s=0 

( j+u+v=4c+m+s ) 


n 

2 


n 

* 


i (d.2_1d'  ® j' g'2G  j ) 
J,k-0  u, m=0  v,s=0  J 11  u ^ s 

( j+u+v=k+m+s ) 


(5-8) 
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Also  note  that  for  fixed  n,  ^wtn^t=p+l  forms  a sequence  of 

m=(i+2n+q)  - dependent  random  vectors,  and  the  W^n  possess 
uniformly  bounded  fourth  moments,  since  the  have  finite  fourth 

order  moments.  Thus  it  follows  by  Theorem  1 in  Schonfeld  (1971) 
that 

1 T 

Zm  = Z W.  has  a limiting  multivariate  normal  distribution 

Tn  yT  t=p+l  tn 

as  T-*oo,  with  mean  vector  0 and  covariance  matrix  V . Then 
using  a vector  analogue  of  Theorem  7-7.1  in  Anderson  (1971),  as  given 
by  Lemma  2 in  Schonfeld  (1971),  we  can  conclude  that 


1 


(lra»^iY)'G,"1(Z“1®  IN)e 


ZTn  + ^ 


converges  in  distribution  as  T->«o  to  a multivariate  normal 
distribution  with  mean  vector  0 and  covariance  matrix  equal  to 


lim  V 
n 00 


00  00 

z z 

J , k=0  u, m=0 


v, s=0  0 


( j+u+v=k+ir.+s ) 


= J G"1Z_1G*“1  ® (A*"1G*ZGA“1) ' dX 

2tt  -tt 

= Yq  . (3-9) 

A similar  argument  can  be  used  to  show  that  all  the  subvectors  of 

— W,G,‘1(Z"1  ® Iv)e 

VT 

have  a joint  limiting  distribution  as  T^-oo  which  is  multivariate 
normal  with  mean  vector  0 and  covariance  matrix  V,  and  hence 
the  theorem  is  established. 
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4.  EXTENSIONS 


The  theory  and  method  of  estimation  in  this  paper  can  easily 
be  extended  to  the  case  of  the  mixed  vector  autoregressive  moving 
average  model  with  exogenous  variables,  the  so  called  VARMAX 
model,  of  the  form 

Yt  - j/iVi  = j0BlXt-i  + €t  + j^t-i’ 

where  the  observable  exogenous  variables  X^  are  suitably  prescribed. 

In  fact,  the  method  which  has  been  used  here  for  the  estimation  of 
the  vector  ARMA  model  has  been  applied  by  Reinsel  (1976b)  to  the  more 
general  dynamic  simultaneous  equations  model  with  autoregressive  moving 
average  disturbances, 

-Ao>'Yt  - l/iVi  - j0Bixt-i  + ut>  Ao  * °- 

where  the  disturbances  satisfy 


Ut  ‘ j^t-i  = €t  +J1Gi€t-i' 
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